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Abstract. The divergence of the constraint quantities is a major problem in computational grav- 
ity today. Apparently, there are two sources for constraint violations. The use of boundary con- 
ditions which are not compatible with the constraint equations inadvertently leads to 'constraint 
violating modes' propagating into the computational domain from the boundary. The other source 
for constraint violation is intrinsic. It is already present in the initial value problem, i.e. even when 
no boundary conditions have to be specified. Its origin is due to the instability of the constraint 
surface in the phase space of initial conditions for the time evolution equations. In this paper, 
we present a technique to study in detail how this instability depends on gauge parameters. We 
demonstrate this for the influence of the choice of the time foliation in context of the Weyl system. 
This system is the essential hyperbolic part in various formulations of the Einstein equations. 



1. Introduction 

One of the major problems in computational gravity is the fact that the constraints are not pre- 
served in free evolution codes. Indeed, it can be observed in many numerical approaches that 
the constraints are violated with an exponential (or even worse) rate in time. Thus, the numer- 
ically generated solution of the evolution equations ceases to satisfy the full Einstein equations 
as time progresses. 

Currently there are two known sources for constraint violation in an initial-boundary-value 
problem. The first one is present already in the Cauchy problem. It is due to the structure of 
the field equations and the specific splitting of these equations into evolution and constraint 
equations. The other cause for constraint divergence is due to inappropriate boundary condi- 
tions, i.e., data given on the boundary of the computational domain which are not compatible 
with the constraint equations. These data will give rise to 'constraint violating modes' which 
propagate into the computational domain thereby spoiling the solution inside. 

Much work has gone in recent years into the possibilities of curing the desease of diver- 
ging constraints. There have been various proposals for constraint preserving boundary condi- 
tions 0|llll||l5l|l6l|28l|3Ol|3lI|32|to prevent the constraint violating modes from entering the 
computational domain. However, the only formulation of an initial-boundary-value problem 
for the Einstein equation which is known to be well-posed has been given by Friedrich and 
Nagy 1 12 1 . On the other hand, if the constraints have already started to diverge, there are ways 
to force the solution back onto the constraint surface I2"l l21ll55lll8l . 

In the present paper we want to discuss the first cause of constraint violation which is re- 
lated to the structure of the field equations. Our aim is to demonstrate that the stability of 
the constraint propagation depends heavily on the choice of coordinates, in particular on the 
time foliation. We will do this on the basis of an explicit example, the Bianchi equation. This 
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equation features prominently in various formulations of the field equations of general relativ- 
ity 1 7 12, 10 1 where it can clearly be seen that it is the essential hyperbolic equation in general 
relativity It is the equation which governs the propagation of the gravitational degrees of free- 
dom described by the Weyl curvature tensor. However, we want to stress that the method we 
employ is general and can be applied to any system of constrained evolution equations. In fact, 
a similar analysis can and should be carried out for the standard Einstein equations in the ADM 
or BSSN formulations. 

The plan of the paper is as follows. In section we present our geometric point of view 
and discuss our approach in more detail. In section [3] and |1] we indicate how to derive the 
evolution and constraint equations for the Weyl tensor and how to find the subsidiary system 
of propagation equations for the constraints. In section|5]we apply the Routh-Hurwitz criterion 
for stability to a suitably simplified set of equations and evaluate the ensuing conditions. Since 
we are not able to give complete mathematical proofs for all the statements made, we dicuss 
at the end of that section the numerical evidence for our claims. We end the paper with a brief 
discussion of the results and implications for further studies. 

2. General description of the method 

In the situation, we are considering, we have to deal with fields on space-time which constitute 
a system with infinite dimensions. In order to get a feeling for the geometrical situation, we 
pretend that we are only concerned with finite dimensions. The material in this section serves 
mainly as a motivation for the calculations performed in the following sections. It is not essen- 
tial for the rest of paper. So let V denote a finite dimensional manifold which we call the phase 
space of the system. We assume that there is a vector field V defined on V , whose integral curves 
describe the evolution in time of the system from some specified initial condition p £ V . Thus, 
V can also be interpreted as the manifold of initial conditions for the system. 

Let O : V — > Z be the constraint map, mapping the initial conditions onto the constraint 
quantities which form a manifold Z. Consider the equation z = ®(p) and assume that there is 
Po £ V and Zq £ Z such that Zq = ®(po)- Let us assume that dCD(p ) is surjective. Then we 
can locally write p = (q, r) so that 3® /dq is invertible. By the implicit function theorem we can 
solve the equation z = 0(q, r) locally near po = (q0/ r o) i- e -/ we find a smooth map q = ip(z, r) 
with qo = i/>(zo,ro) such that z = ®(t/>(z, r),r) for all (z, r) close to (zq, tq). This allows us to 
locally consider the phase space V as being parameterised by p = (xp(z, r), r) where z are the 
constraint quantities and r comprises the 'residual' variables. These are the unconstrained 'true 
degrees of freedom'. Thus, V is locally foliated by the leaves of constant z. 

Consider now the vector field V on V . Given an arbitrary point po 6 V there is a unique 
integral curve pt of V through po such that pt = V(pt). Related to this curve is the curve 
Zf = O(pt) in Z. This curve describes the change in the constraint variables caused by the 
evolution. Its tangent vector it = dO(pf) ■ pt = d<$)(pt) ■ V(pt) depends on the solution curve 
pt- Using the parameterisation for pt we can write 

(2.1) z t = F(zt,rt), 

for some smooth map F : V — > TZ. Now it may happen that for some zq £ Z we have 
F(zQ,r) = for all r. Then zq = O(pt) for all times t if zq = O(po), i.e., the evolution remains 
in the leaf C = O~ 1 (zo), the constraint manifold. In that case, one says that the constraints are 
propagated by the evolution. This is the case for the Einstein system in its various formulations 
and for many other constrained evolution systems appearing in physics (however, see 1 8 [ for a 
well-known example of a system where this is not the case). 

We are interested in the case when the solution curve starts outside of, but close to, C and we 
want to obtain some information about the change of z during the evolution. So let us take a 
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1-parameter family of evolutions pt (A) with the corresponding 1-parameter family of constraint 
variables Zf(A). Assume that Zf(0) = Zq so that Pf(0) = (zq, ff) lies on the constraint surface. 
Then we have 

(2.2) ^(A)=F(z f (A),r t (A)). 
Taking the derivative with respect to A and evaluating at A = yields 

(2.3) tt-<5z = F z (z ,r t ) ■ 5z + F r (z ,r t ) ■ 5r = F z (z ,r f ) ■ 5z. 

The second equality follows by taking the derivative of F(zq, r) = with respect to r. This 
equation tells us how perturbations in the constraint variables close to the constraint surface 
evolve once they are excited. Their propagation properties are determined by the evolution pf. 

In the context of the Einstein equations the evolution curve pt corresponds to the space-time 
which evolves from the initial conditions po- Thus, the perturbations of the constraint quantities 
propagate on the background space-time provided by the solution under consideration. 

In the finite dimensional setting this is a rather straightforward route to determine the propaga- 
tion properties of constraint perturbations near C. However, the Einstein system is infinite 
dimensional and it is not so clear how much of this route translates rigorously to an infinite- 
dimensional setting. In a recent paper, Bartnik |1 J shows that much of the finite-dimensional 
picture can be taken over to the Einstein system on asymptotically flat manifolds in the for- 
mulation given by Fischer and Marsden 1 6 1 . In particular, he shows that the constraint map 
O is smooth and surjective and that all its level sets, in particular the constraint manifold C, 
are smooth Hilbert submanif olds of the phase space of GR defined by the first and second fun- 
damental forms (g, n) of a suitable 3-dimensional manifold. Thus, the infinite dimensional 
Einstein system shows some of the features as the finite-dimensional model described above. 

The leap from the finite-dimensional to the infinite-dimensional system is quite considerable. 
In order to make contact with what follows let us remark that the map F will in general contain 
partial derivatives of the constraint quantities. It is a general procedure to reduce PDEs to 
ODEs by Fourier transform. However, in order to apply this in the present situation several 
assumptions need to be made (see app.|A]). Taking the finite-dimensional case as motivation, 
we are therefore led to study the linearisation of the system which propagates the constraints. 
For a lack of a better name we call this system following Friedrich the subsidiary system. Since 
this system is linear in the constraint quantities we have to study this system itself. In the 
course of the investigation we may assume that the background manifold is a fixed solution of 
the Einstein equations. Of course, this procedure is not limited to this particular formulation of 
the Einstein equations but applies to any formulation for which the constraints are propagated 
by the evolution equations. 

In fact, following Frittelli |13[, Shinkai and Yoneda |24||55J 157||33 have already studied 
the stability properties of the subsidiary system for several variations of the ADM system, most 
notably the BSSN formulation. Their work has been motivated by the desire to understand 
the superiority of the BSSN scheme over the standard ADM formulation. They analyse several 
modified ADM formulations on flat space or on the Schwarzschild space-time with a fixed time 
foliation. Compared with their approach, our work will be both more restrictive and more gen- 
eral. We do not restrict ourselves to a given background and admit arbitrary time foliations. But 
since our aim is to determine the propagation properties analytically, we cannot easily switch 
between various different formulations because the algebra is rather complicated. 

In this paper, we will start with the analysis of the system of constraint propagation for a 
class of formulations of the field equations given by Friedrich 191 1101 1111 following this route. 
These are first order formulations consisting of equations for several geometrical quantities, 
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most notably the extrinsic curvature and the acceleration vector of the time foliation and various 
curvature quantities including the Weyl curvature. 

The size and resulting complexity of the system makes it hopeless to analyse the full system 
at once. Therefore, we restrict ourselves to the most important subsystem which features in 
all these formulations, namely the so-called Bianchi equation. This is an equation for the Weyl 
curvature tensor which results from the Bianchi identity for the Riemann tensor after using the 
Einstein vacuum equations. In the conformal setting there is an additional conformal rescaling 
involved which, however, does not change the character of the equation. 

As will be described in more detail in the following sections, the Bianchi equation can be 
split in the usual way into constraint and evolution equations, and it can be verified that the 
constraints propagate. This collection of constraints and evolution equations will be referred to 
as the Weyl system. The subsidiary system of evolution equations for the constraint quantities 
contains not only the constraint quantities of that system but also constraint quantities whose 
vanishing implies the consistency of the acceleration and extrinsic curvature with the existence 
of a foliation. These constraints arise when acceleration and extrinsic curvature are also evolved 
numerically. We decouple the Weyl subsystem from the other subsystems by assuming these 
constraints to vanish identically i.e., that only the perturbations of the Weyl constraints are 
excited. This amounts to the assumption that the Weyl curvature propagates on a foliation 
which is not influenced by the curvature and vice versa. 

Having obtained the subsidiary system which is already linear in the constraint quantities, 
the next step is to localise the equation by 'freezing' the coefficients. This means that we study 
the system in an infinitesimal neighbourhood of an arbitrary but fixed event. This results in 
a system with constant coefficients which can be treated by Fourier analysis. We derive the 
mode dependent propagation matrix P(k) and ask for its stability properties. The main tool 
in this analysis is the Routh-Hurwitz criterion which allows us to determine the number of 
eigenvalues of P(k) with negative real part by looking at the coefficients of its characteristic 
polynomial. 

One might question the relevance of the frozen problems to the problem with variable coef- 
ficients. This is not an easy task to sort out. One possibility is to refer to the literature on the 
analysis of PDEs such as 1 20 [ where it is shown that for strongly hyperbolic or second order 
parabolic systems well-posedness of all frozen systems is sufficient for well-posedness of the 
general problem provided there exists a smooth symmetrizer. For first-order systems Strang 1 29 1 
has shown that it is also necessary. This indicates that the properties of the frozen systems and 
in particular the estimates which relate the solution at time t to the initial data are closely related 
to those of the general system. 



3. Hyperbolic reduction 

The formulation of an initial-value problem for the Einstein equations, which is the basis for 
their numerical treatment, requires the introduction of a time-flow along which the integration 
of the field equations proceeds to produce a solution out of initial data. The covariant fields are 
decomposed into parts tangential and transversal to this flow ((3+l)-decomposition) which splits 
the originally covariant field equations into a set of equations for the (3+l)-constituents of the 
fields, hopefully yielding a symmetric-hyperbolic system of evolution equations which allows 
for the formulation of a well-posed initial value problem. 

In this section we want to present this procedure and the formalism used here on the example 
of the Bianchi equation: 



(3.1) 



V«K« fecd = 
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The tensor ^ a bcd is a trace-free tensor with the symmetry properties of the Weyl tensor describ- 
ing the gravitational field. Importance and origin of this equation are discussed in 1 23 [. 

3.1. (3+l)-decomposition. We work with the time-like, normalised vector field t" generating 
the time-flow and use metric signature (+,—,—,—), thus t"t a = 1. With respect to this vector 
field, every tangent space splits into a parallel (1-dim. time-like) component and an orthogonal 
(3-dim. space-like) component. The respective projectors are 

(3.2) t a b = t% and h a b = 6 a b - t% , 

which also splits the metric: 

(3-3) gab = tah + h a b, 

where h ab is the negative definite spatial metric in the space transversal to t a . 

Accordingly, every tensor splits into parts which are parallel or transversal to t a . We call 
a tensor purely spatial if every contraction with t" or t a vanishes. As an example, the (3+1)- 
decomposition of K a bcd is 

(3-4) K abcd = 4t[ a E b ^ c t d ] + 2tyB b ^t ecd - 2£ ab B e ^ c t d ^ + e abe E e Uf cd 

with the purely spatial, trace-free and symmetric tensors E ab and B ab which are called electric 
and magnetic component of the gravitational field. The covariant derivative is decomposed as well: 



(3.5) V„ = t a D + D a 
with its components 

(3.6) D := t a V a and D„ := h h a V b . 

To facilitate calculations, it is useful to introduce derivative operators which are adapted to the 
time vector-field t" 1 7 p. 65] . To characterise the course of t a , we introduce the extrinsic curvature 
quantities 

(3.7) X a -=Dt a and X a ■= D a t b (with trace X := X a ") ■ 

Since t a is assumed to be normalised, these quantities are purely spatial. Usually t a is chosen 
as the unit normal field of a foliation of space-time. Then Xab is the extrinsic curvature of the 
foliation and symmetric. Taking t a as the 4-velocity of an observer travelling along the integral 
curves of t" , then x" = t b X7t,t" is the acceleration measured by the observer. Therefore x" is 
called acceleration vector. The adapted derivatives are defined by their action on 1-forms: 

(3.8) dv b := Dv b + t b x a v a - Xbt a v a , 

(3.9) d a v b := D a v b + t b x a c v c - Xab t c v c 

Their action on vector fields and higher tensors is defined by the Leibniz rule and the require- 
ment that when applied to functions, they coincide with D and D a . 

The adapted derivatives have the important property, that in contrast to D and D a they com- 
mute with the projectors defined in <3.2t . The new time-derivative 3 can further be interpreted 
as the generator of Fermi- Walker transport along the integral curves of t a . The spatial derivative 
d a is the Levi-Civita connection intrinsic to the leaves of the foliation. 
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3.2. The Weyl system. With the mentioned tools, we are in a position to carry out the (3+1)- 
decomposition of the Bianchi equation 

ED V a K a bcd = 0. 

We first decompose the covariant derivative according to J3.5I , then transform to the new de- 
rivatives d and d a by use of J3.8I and J3.9I and finally insert the (3+l)-representation of the 
gravitational field as given by 13.41 . The resulting equation still has three indices but each of its 
terms can now easily be classified to be either purely temporal or purely spatial in any of its in- 
dices. This requires every such component of the equation to hold on its own, thereby splitting 
the equation into a set of equations for the (3+l)-components E ab and B ab of the gravitational 
field. 

The calculations are too technical and lengthy to be given here 1 34 1, therefore we only give 
the resulting Weyl system of equations: 
The constraint equation for E: 

(3.10) d a E ac = e ahe Xab B ec + e cbeXa h B ae 
The constraint equation for B: 

(3.11) d a B ac = -e abe Xab E ec - e cbeXa h E ae 
The evolution equation for E: 

(3.12) dE bc + e ae{h d a B c) e = -2 X E bc + 2 X \ b E c)a + X[c a E b)a - h bcX ae E ae + 2 X a e ae[b B c) e 
The evolution equation for £>: 

(3.13) dB bc - e ae{b d a E c) e = -2 X B bc + 2 X \ h B c)a + X(c a B b)a - h bcX ae B ae - 2 X a e ae{b E c) e 

The system of equations has the remarkable property that it is invariant under the duality trans- 
formation E — > B, B — > — E known from electrodynamics. This allows for a very compact and 
elegant notation: Formally collecting E and B into a complex tensor F bc := E bc — iB bc , the dual- 
ity transformation now becomes simply F bc — ► i F bc . With this, the set of equations reduces to 
one single (complex) constraint equation 

(3.14) d a F ac = ie abe Xab F ec + U cbeXa h F ae 
and one single (complex) evolution equation 

(3-15) dF bc + ie ae{b d a F c) e = -2 X F bc + 2 X a (b F c)a + X{c a F b)a 

-h bcX ae F ae + 2i X a e ae[b F c) e . 

It can be shown 1 9 1 that the evolution equation is symmetric-hyperbolic and that therefore its 
initial value problem is well-posed. The development of the gravitational field is completely 
determined by the evolution equation. Thus it defines the vector field V in the picture of sect.|2] 

4. Constraint propagation 

The last statement gives rise to the question of compatibility between the evolution and con- 
straint equations: Given data on an initial time slice which fulfill the constraint equation, then 
the evolution equation fully determines the time development of these data, but will the con- 
straint equation hold on later time slices as well? In other words: Is the evolution-generating 
vector field V tangential to the constraint surface C? 

This important question can be answered by looking at the time development of the con- 
straints. Therefore, we write the constraint equations (13.101 and <I3. 1 II as 

(4.1) = £ c := d a E ac - e abe Xab B ec - e cbeXa h B ae , 
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(4.2) = B c := d a B ac + e ahe Xab^ec + e cbe x a b E ae 

with the constraint quantities £ c and B c whose evolution equations we need to determine. Due 
to the invariance of the system of equations under duality transformation we need to calculate 
only one of them. The other one is then obtained by substituting E — > B, B — > — E and £ —> B, 
B-+-E. 

According to <4.1> . calculating d£ c gives on the right hand side time derivatives of B for 
which its evolution equation can be substituted, but furthermore the time derivative of a spatial 
divergence of E. Using the evolution equation of E in this place makes it necessary to commute 
the spatial and time derivative producing curvature terms. Specialising to the case that t a is in 
fact orthogonal to a foliation and thus assuming Xab to be symmetric, finally yields the following 
system of propagation equations: 

(4.3) d£ c = -\e cah d"B h + \e cahX a B h - \ X £ C + & c % 

- E ab K-cab + 2E ac K. ha b + 2e ab ( c B d b lC da , 

(4.4) dB c = \e cah d a £ b - \e cahX a £ b ~ \xB c + \x?®b 

- B ah K, cab + 2B ac )C ba b - 2e ab ( c E d b K, da 

with the constraint quantities of the foliation: 

(4.5) K ab := 3[ fl x fc] , 

(4.6) K, abc := 23lV ]c + h a a ,h w h cc 't d 'R a ' Vc , d , 

Vanishing of JC ab guarantees Xab to remain symmetric during evolution, vanishing of K, abc 
is equivalent to the second Gauss-Codazzi relation which has to hold if Xab is the extrinsic 
curvature of a foliation. 

The first lines in J4.31 and J4.4> feed back the constraint quantities £ and B into themselves 
with the extrinsic curvature quantities of the foliation acting as coefficients. The second lines 
couple the system to the constraint quantities of the foliation with E and B acting as coefficients. 

Obviously all the terms on the right-hand sides are proportional to constraint quantities. 
The differential equations therefore are homogeneous with respect to constraint quantities, i.e. 
£ c = 0, B c = are solutions of the propagation equations under the assumption, that the 
constraint quantities of the foliation also vanish. Since the system can be shown to be symmetric 
hyperbolic, we have uniqueness of solutions, so that the only solution with vanishing initial 
conditions is in fact the zero solution. Hence it is shown that the evolution and constraint 
equations of the Weyl system are compatible. 



5. Stability analysis 

From the analytical point of view, the above statement is all we need. From the numerical point 
of view, this is just the first step. 

In doing numerics, the canonical approach is to solve the constraint equations to produce 
initial data which is as accurate as possible. The numerical integration of the evolution equations 
will produce a solution of the evolution equations from these initial conditions. In general, this 
procedure cannot distinguish between good data satisfying the constraint equations exactly and 
bad data which are perturbed by numerical error. Numerical noise will be carried along and 
can (and will) accumulate. 
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That means that one has to consider the evolution equations from a more general point of 
view allowing arbitrary initial data (off the constraint surface C). The way in which these non- 
solutions of the constraint equations are propagated, depends on properties which are not part 
of the original full system of equations but of the evolution equations on their own. 

The form of this system is not fixed. By adding multiples of the constraint equations one 
can write down many different systems with (presumably) very different properties, changing 
the evolution anywhere but on the constraint surface C . Here we will consider the form of the 
evolution equations as fixed, partly because we want to focus on the particular influence of the 
foliation and partly because the spinorial formulation of the Weyl system seems to suggest that 
this form is a very natural one. 

Then we see that the subsidiary system J4.3I4.4I essentially depends on the foliation chosen 
for hyperbolic reduction as this is the parameter which determines how the properties of the 
covariant equation <3.1> are partitioned between evolution and constraint equations. 

Since the numerical procedure aims at producing solutions which are as close to analytic 
solutions as possible, it must be required to be stable against perturbations by numerical noise. 
Thus, it is a necessary condition that the solutions of the evolution equations with vanishing 
constraint quantities are attractors in the positive time direction, i.e., the constraint surface in 
the phase space of initial conditions has to be attractive. 

The following analysis will extend the analysis of compatibility given in the last section, 
which can be looked upon as stability analysis of zeroth order, to a stability analysis of first 
order which is valid for small perturbations. In this process different approximations have to 
be applied. The first one is that we analyse the constraint propagation properties only within 
the Weyl system. The coupling to other equations outside this system will be neglected by 
assuming the external constraint quantities (belonging to equations for the foliation) to vanish. 

To make the calculations more compact, we introduce a complex constraint quantity T c := 
£ c — iB c to exploit the invariance under duality transformation. Then the decoupled propagation 
equations combine into a single one which reads: 

(5.1) dT c = d£ c - idB c 

_ 1 ; f a« -r-b _i_ 3.V v fl fb i 1 h _3 v -r- 

1 lt cab IJ •> ~r 2 L cabX •> < •'b 2"- c 

Obviously, T c = is a solution, but now of particular interest is, how solutions T c 7^ will 
behave. If propagation is stable, they will converge against T c = in positive time direction. If 
not, then the constraint quantity will diverge. 

To investigate this, we use another approximation: In general, the coefficients of the propaga- 
tion equation Xa an d x" vary from point to point. We now consider the propagation properties 
locally around a point p 6 Ai and assume, that in a certain neighbourhood of this point the 
coefficients can be considered constant. That implies, that the space-time manifold is locally 
approximated by its tangent space at point p. This will be the manifold of our further investig- 
ation. The constraint quantity T c and the extrinsic curvature quantities Xa an d X a accordingly 
become tensor fields on the flat tangent space, for which now is imposed a differential equation 
with constant coefficients which formally corresponds to the original equation. The detailed 
discussion of what is involved in this step is given in the appendix lAl There we show how to 
derive the final equation <A.5t which will be analysed here. Note, that this equation contains 
the lapse function N and the shift vector. However, for the purpose here, it is enough to assume 
N — 1 and a vanishing shift, so that d = df. We will comment on the influence of non- trivial 
lapse and shift below. 

The procedure of this approximation is known as freezing of coefficients, a standard method 
in numerical stability analysis. Because the frozen propagation equation is now defined on flat 
space, it is apt to be Fourier transformed in the spatial directions. Let V denote the local spatial 
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tangent space (tangent to the space-like slice through p) and V* its dual space. Then the frozen 
constraint quantity can be represented as 

(5.2) T c (t,x a ) ■- j f c (t,k a )e ik » xa d 3 k a Vx a eV 

v* 

with its frequency components T c (t,ka), for which now the following propagation equation holds: 

(5.3) d? c = \e cab k a P + - \x?c + \xc?t> 

= V c h {k a , Xa h ,X a )h{t,ka) 

The spatial derivative d a has transformed into the frequency covector k a which reduces the 
propagation equation to an ordinary differential equation with constant coefficients which are 
denoted by the propagation tensor 

(5.4) V c h := e ca b (\k a + \i X a ) + \ Xc b - \ X h h c . 

The propagation tensor is the evolution generator of the constraint quantities in frequency rep- 
resentation and its eigenvalues decide on the propagation properties of the mode belonging to 
the respective eigenvalue. 

For numerically stable constraint propagation, the propagation equation is required to be 
stable and attractive around the point T c = 0. Here stable means that the solutions around 
J- c = can be controlled: For every maximal deviation from J- c = given for all times, there 
is a maximum initial deviation. Attractive means that in a certain neighbourhood of T c = 
all solutions converge against T c = for large times. If both conditions are met, then the 
propagation equation is said to be asymptotically stable which is equivalent to all eigenvalues 
lying in the left complex half -plane {£H(z) < 0}. More details can be found in the literature on 
linear systems, e.g. |5|. The impact of the location of the eigenvalues and of diagonalisability 
on the propagation properties has been analysed by Yoneda and Shinkai in |35|. 

To study the eigenvalues, it is necessary to calculate the characteristic polynomial 

(5.5) X v (z) := det(V a h - zh b a ) , 

which is most elegantly done by using the covariant representation of 3-dimensional determin- 
ants as 

(5.6) detT fl b = T fl [ fl T fc fc T c c ]. 

Since the propagation tensor is only 3-dimensional in our case, it would be possible in principle 
to directly calculate the eigenvalues by the well known solution formula for the roots of poly- 
nomials of third order. Unfortunately, the solution formula employs case discriminations which 
makes the general dependence of the eigenvalues on the parameters of the propagation tensor 
difficult to analyse. Moreover we are only interested in the sign of the real part of the eigenval- 
ues to decide whether propagation is stable or not. Thus, we only need to know under what 
conditions the spectrum of the propagation tensor is contained in the left half of the complex 
plane. This is exactly the kind of question that can be decided with the Routh-Hurwitz criterion 
(see app.|Bj which is applicable to propagation tensors of arbitrary dimensionality. 

5.1. Application to the propagation tensor. The first step now is to calculate the characteristic 
polynomial in the representation required by the Routh-Hurwitz criterion and results in 

(5.7) X-p{iz) = b z 3 + biz 2 + b 2 z + b 3 +i (a z 3 + a x z 2 + a 2 z + a 3 ) 
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\{kak a -9 X aX a ] 



with the following real coefficients: 
(5.8) 8 = 1, «i = 0, 

a 2 = -(\{M a a f-\M h c M c 

a 3 = -\k a M a h Xhl 
b = 0, b 1 = -M a a , 
b 2 = -^k aX a , 

b 3 = detM a b - \M a b (k% - 9 X a Xb) 

Here, M a \, denotes the symmetric part of the propagation tensor which is a trace-transform of 
the extrinsic curvature: 



(5.9) 



M 



cb ■- 



2 Xcb 



cb 



From these coefficients we calculate the three Hurwitz determinants T> 2 , T> 4 and T>(, which are 
required to be strictly positive for stable constraint propagation. The three inequalities are of 
increasing complexity since 23; is a i x /-determinant. 

5.1.1. The first Routh-Hurwitz inequality. The condition requiring T> 2 to be strictly positive, is 
(5.10) < V 2 = 





«i 




1 





bo 


h 







-M/ 



-M " = 4y 

a m 



(RH.l) 



& X>0 



Remarkably, this condition neither depends on the mode k a , nor on the acceleration vector x" ■ 
It demands that the foliation has positive mean curvature at the point under consideration. 

5.1.2. The second Routh-Hurwitz inequality. This inequality is already much more complicated, 
therefore we have to limit ourselves to the results. The complete calculations can be found in 
1341 sec. 6.3]. The determinant involved here is: 



(5.11) 



X>4 



fl «i a 2 a 3 
bo h h b 3 

flg fll a 2 







(\k aX a ) 



bo b\ b 2 

= \M C C (detAC fc - \M a h (k% - 9 X a Xb), 
with another trace-transform of the extrinsic curvature 

(5.12) M b := M a h - M c c h h a = \xa + h h a > which implies 

15.41 

(5.13) M a a = M a a -3M/ = -2Mf = 8 X . 



This determinant now contains all the parameters k a , x a and Xa ■ Ideally one would like to fulfil 
this condition for arbitrary modes k a simultaneously, resulting in relations between x" an d Xa 
only. Detailed analysis shows that this is actually possible and yields the following conditions: 
Let the acceleration vector x" be represented in polar fashion as x" = s v " with an unit vector 
v a and s > 0. Then the following conditions are necessary and sufficient for the second Routh- 
Hurwitz inequality D4 > to hold for arbitrary modes k a : 

(RH.2) All three eigenvalues n,- of AfJ 1 are strictly positive; 
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for the length s of the acceleration vector hold both 

2 



(RH.3) 



(RH.4) 



s < 



s < 



detfif a b 

M ah v a v b 



and 



Here Af a b denotes the inverse of Af a b . The condition JRH.2> implies the first stability condition 
JRH.H . The conditions JRH.3> and JRH.4> are not equivalent. Examples show, that depending 
on the parameters, either one of them can be more strict than the other. 

5.1.3. The third Routh-Hnrwitz inequality. The third inequality is obtained from the 6x6- 
determinant 



(5.14) 



-D 6 



ciq a\ a 2 a 3 

b bi b 2 b 3 

«o a l a 2 a 3 

b b x b 2 b 3 

«o a i a 2 





3„ 2 









«3 

b\ b 2 b 3 
2„ 2 



a 3 — bi a 2 b 3 + b 1 a 2 b 2 a 3 + 2b\a 2 b 3 — 3bib 2 a 3 b 3 

- a 2 b 2 2 b 3 + b 2 3 a 3 - b 3 , 

with the coefficients already given in Il5.8> . To examine T>^ > with respect to the frequency k a , 
we sort the terms in by their order in k a : 

(5.15) V 6 =K 6 + K 4 + K 2 + Kq 
The terms X, contain k a to the power of i and are given by 

(5.16) K 6 :=B 1 ((B 1 -fo 1 A 1 ) 2 -fc 2 2 A 1 ), 
(B 1 - b 1 A 1 )(2b 1 B 1 {A + A 2 ) + (B + B 2 ){A 1 b 1 - 3B t ) + b 1 b 2 a 3 ) 



(5.17) 
(5.18) 



Hi 
X 4 



+ b 2 2 Bi(A + A 2 ) + b 2 2 Ai{B + B 2 ) + lb 2 a 3 B x b x + b 2 3 a 3 , 



b 1 2 B 1 {A + A 2 ) 2 
- (A + A 2 )(B 



h(B + B 2 ) 2 (3B 1 -2A 1 b 1 ) 
B 2 )(bi(4Bi -2h A x ) + b 2 ) 



+ b 2 a 3 b x 2 (A + A 2 ) - 3b 2 a 3 b 1 (B + B 2 ) - b x 3 a 3 2 , 



(5.19) K := -(B + B 2 )((B + B 2 )-b 1 (A + A 2 ))\ 
with the following abbreviations: 

A := -M} a M b h KA 1 := -\k a k a ,A 2 := -\x a Xa) 

B := := \M*k a k h ,B 2 := \M^Xb 

This makes it obvious that we have to discuss a polynomial inequality of sixth order in k a of 
which we hope, that it can be fulfilled for arbitrary k a simultaneously with the first and second 
Routh-Hurwitz inequalities above. 

To see, if this is actually possible, we will first discuss the low frequency limit. In the domain 
of low frequencies, the terms of low order in k a play the dominant role. Looking at the case 
k a — > 0, only — > K$ contributes. Analysing this term yields the following result: 
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As before let A4 a h denote the symmetric part of the propagation tensor and choose the ac- 
celeration vector x" to be represented in polar form as x" = s v " with unit vector v" and length 
s > 0. 

Then necessary and sufficient for the third Routh-Hurwitz inequality T>(, > to hold in the 
limit k a — > 0, are the following conditions: 

(RH.5) All three eigenvalues m, of Ai^ are strictly negative ; 



(RH.6) 




det.M fl b 
M a i,v a v b ' 



Because of = m.{ — {m\ + ni2 + mo l ) = — Y,j^i m j, it follows from JRH.5> . that n, > ^RH.2> 
which further implies (|RH.1> . Moreover, it can be shown, that <RH.6> implies the previous 
conditions <RH.3> and (|RH.4> . Therefore the conditions JRH.51 and <RH.6> are already sufficient 
for the first and second Routh-Hurwitz inequality T>2, £>4 > to hold, and that for all k a . 

For non-zero frequencies, the inequality X>g > 0, of course, is much more complicated to 
analyse. Representing the frequency vector in polar form as k" = k u" with another unit vector 
u a and k > results in: 

(5.20) V 6 = K 6 {u a )k 6 + K 4 {u a )k i + K 2 {u a )k 2 + K {u a ) 

Surprisingly, numerical tests (see below) strongly support the conjecture, that K(,(u a ), K^(u a ) 
and K2(u a ) are already individually positive whenever the conditions JRH.51 and <RH.6> hold. 
This means, that zero-frequency stability already implies stability for arbitrary modes. 

To test this conjecture in a systematic way, it is useful to represent the acceleration vector x" 
in the following rescaled polar form (with gam t > and direction v a v" = —1): 

(5.21) x"(£,u fl ) := tX a (v a ) with 

(5.22) X>«) := I 



I / detM* ^ 

M a bV a V b 

Then JRH.61 is equivalent with f < 1. Inserting this into the coefficients Xg, K4 and K2 of 15.201 
has the following benefit: Rescaling M. b — > aM. b now causes K, to rescale with a power of a: 
Ki -> a 9 - ! 'X/. That means, that rescaling with positive factors will never change the sign of the 
individual terms. 

Without loss of generality assume that the eigenvalues m, of M. b are numbered in increasing 
order: m\ < < < 0. Then any Ai a b can be represented as Ai a b = —m\ Ai a h , where the 
eigenvalues of M. b are given by —1 = fh\ < ih^ < m$ < and the X,- will have the same 
signs for both tensors with and without tilde. This shows, that it is sufficient to prove or test the 
conjecture for eigenvalues ( — l,ni2,nii,) lying inside the bounded triangle defined by the last 
inequality. 

Although we did not find a proof for the given conjecture, we neither found any special cases 
in numerical tests, in which the conjecture would be falsified. For the numerical tests, we used 
the following procedure: First, pick the following quantities: 

• m.2, mo, with — 1 < m.2 < M3 < 0, 

• the gain t with t < 1, 

• the direction v a with v a v a = —1, 

• the direction u a with u a ii a = —1. 

These amount to seven real and independent degrees of freedom, all bounded to finite ranges. 
The conditions of the first two items assert, that the conditions JRH.5> and JRH.6> hold. 
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Second, calculate the coefficients K^, K± and K2 of ( 15.201 . If they are all positive, then the 
conjecture is proven for the chosen special case. If K(, is negative, the conjecture is falsified for 
the high-frequency limit. If X4 or K2 are negative, a more detailed analysis must show, if this 
results in frequency bands which are instable. 

This procedure has been carried out for a large number of special cases. The eigenvalues 
have been chosen as :— —cc and :— —a/3 with a and /3 in [10~ 5 , 1]. In different runs we 
tried both equally distributed and exponentially distributed values, which concentrate around 
the critical points a = and ft = 0. The typical (a, /3)-grid had 20 x 20 points. 

For the unit vectors v a and u a we chose vectors on the unit sphere in a way, so that the angular 
distance between them stays approximately constant for all latitudes, which is not the case for 
a simple uniform (G, $)-grid. So the number of points on each line of latitude increases towards 
the equator. Further it was taken advantage of the fact, that the quantities K( are invariant under 
inversion of the unit vectors, so only the northern hemisphere was actually used. On typical 
runs, we used 272 points on this hemisphere (10 latitudes, 40 longitudes on the equator). 

For the choice of the gain t, which describes if the condition <RH.6> is met, we also used an 
exponential spacing, which concentrates around the critical value t — 1. We chose 20 values in 
[0, 1 - 10" 5 ]. This amounts to 591.872.000 combinations per run. 

The result of this test is, that in the given domain, we did not find any cases contradicting 
the conjecture. Only when choosing values t > 1 or a, j3 < we found negative coefficients Kj. 
As the coefficients Kj only consist of (admittedly complicated) polynomials, we tend to believe 
that the conjecture is true. 



6. Results 

In the last section it has been shown that it is possible to apply the Routh-Hurwitz criterion to 
the propagation tensor and to perform a detailed analysis of the geometrical meaning of the 
resulting conditions. In our example of the Bianchi equation, this results in: 

For the extrinsic curvature Xab define the auxiliary tensor M. a \, :— \xab ~ \xf n ab- Represent 
the acceleration vector (see 13.11 in polar fashion as x" — s v " with unit vector v a and positive 
length s. Then the following conditions are necessary for local stability: 

RH.5: All eigenvalues m; of Ai a b are strictly negative; 

RH.6: the length s of x" is, depending on its direction v a , bounded by 



s < 



2 / detM a h 

3 V M ab Vv b ' 



Here and in the following local stability means asymptotic stability of the problem which results 
from localisation in the sense of the freezing of coefficients approximation. 

According to the discussion at the end of 15.1.31 we conjecture, that these conditions are also 
sufficient for local stability. It should be stressed, that these specific conditions of course only 
apply to numerical calculations which explicitly integrate the Weyl subsystem analysed here. 

The conditions in a way resemble the partition of the original set of equations into constraint 
and evolution equations: As presented in 13.11 the extrinsic curvature Xab can be defined as the 
purely-spatial derivative of the foliation's normal unit vector field t a , whereas the acceleration 
vector x" is its temporal derivative. As JRH.5> only contains the extrinsic curvature, it can be 
viewed as a constraint inequality, because it has to hold on every individual leaf of the foliation. 
On the other hand, <RH.6> connects the extrinsic curvature with the acceleration vector, thus 
forming an evolution inequality for the vector field. 

We have derived these conditions under the assumption of a constant lapse and a vanish- 
ing shift. It is straightforward to incorporate the case of non-trivial lapse and shift as follows. 
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Consider the propagation tensor V b for the general equation <A,5> . It is easy to verify that 

V a b = NP a b + i(N l k l )h a b . 

Thus, the effect of the lapse on the spectrum of V b is simply a scaling with a positive number 
while the shift vector shifts the spectrum along the imaginary axis. None of these modifica- 
tions affects the number of roots in the left half of the complex plane. Therefore, the stability 
conditions for the general propagation tensor V b are the same as for V b . 

One may wonder, why there is no influence of the lapse on the propagation properties. After 
all, it is the lapse function which determines the time-foliation. However, this is easily explained 
since it is not the value of the lapse which is relevant but its spatial derivative and this is related 
to the acceleration vector by 

Xa N ■ 

Thus, the spatial variation of N does have a strong influence on the stability properties. 

Results described by Husa 1 19 1 for Minkowski evolutions with the conformal vacuum field 
equations support the results of our stability analysis. The conformal vacuum equations given 
by Friedrich 1 9 1 contain the Bianchi equation explicitly, and when choosing a static hyperbol- 
oidal gauge which satisfies the gauge conditions in the computational domain, it is observed 
there that even though the metric components show exponential divergence from the analytical 
solution 1 19 figure 7], the curvature invariants / and / which consist of components of the Weyl 
curvature, show exponential decay towards the analytical solution |19 figure 6]. 

6.1. Geometric interpretation. The conditions JRH.l^ . (|RH.2^ and <RH.5> all limit the allowed 
eigenvalues of the extrinsic curvature. Their impact is visualised in fig. ^ Each axis in these 
pictures corresponds to one eigenvalue of the extrinsic curvature. Then only such combinations 
are allowed which lie on the same side of all the depicted planes as the corresponding arrow. 
These pictures show that the conditions are of increasing strictness. 

The strongest condition found for the acceleration vector x a tRH.6l can be rewritten as 

(6.1) (X 1 ) 2 , (X 2 ) 2 , (X 3 ) 2 ^ 4 

f«2ffZ3 miH?3 m\m2 9 

with the components of x" taken with respect to the normalised eigenvector basis of M. b ■ This 
obviously means, that the vector x" has to lie strictly inside an ellipsoid whose semi-major axes 
are determined by the eigenvalues nij. 

6.2. Trivial foliation of Minkowski space. As an example consider Minkowski space with 
standard coordinates (t,x,y,z). Its flat foliation by {t = const} -surfaces with vanishing curvature 
quantities Xab — and x" = does not fulfil the stability conditions: 

Vanishing extrinsic curvature results in M. a \, = which violates JRH5> . However, this viol- 
ation is only marginal: The situation lies exactly on the border of the condition. 

Furthermore, the right hand side of <RH.6> is undefined. The formal singularity can be re- 
solved by taking an appropriate limiting procedure (the numerator is of third order in A4, the 
denominator is only of first order), but this will still only result in the condition s < which 
cannot be satisfied. But also here, the condition <RH6> is only marginally violated. 

So it is to be expected, that the eigenvalues of the propagation tensor lie on the imaginary 
axis. The propagation equation with vanishing curvature quantities reads (15. II 

(6.2) dF c = ~ie cab d a F b , 

which is completely analogous to the Ampere-Faraday law of vacuum electrodynamics. The 
propagation tensor in frequency representation therefore is V b = je ca b k a with its eigenvalues 
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{0, i^^—k a k"}. As expected, their real parts vanish. Therefore flat Minkowski foliations are 
not stable in the strict sense used above, but instead are marginally stable. However, this might 
be unstable enough to spoil numerical simulations, since constraint violations will evolve un- 
damped and so can pile up, even though the growth rate will not be significant. 



a 















FIGURE 1 . Limitation of the eigenvalues of extrinsic curvature by the condi- 
tions a: JRHJ), b: fRH2\ and c: JRRSl 
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6.3. Hyperboloidal foliation of Minkowski space. Now consider the foliation of Minkowski 
space given by {t = const} surfaces of the parametrisation 

sinr sinp 
(6.3) t = , r 



COS T + COS p COS T + COS p 



in spherical coordinates (t, r, 9, (p) with compactified time coordinate |t| < n and compactified 
radial coordinate < p < n — \r\ (see fig.0- The leaves in standard coordinates are shown in 
fig-El 



r 



FIGURE 2. Hyperboloidal foliation of Minkowski space (compactified coordinates) 




FIGURE 3 . Hyperboloidal foliation of Minkowski space (standard coordinates), 
equally spaced in the compactified time coordinate 
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This gives for the extrinsic curvature: 



(6.4) 



(6.5) 



Xa 



V COt 2 T + 1 

4 sign(r) j 



Vcot 2 T + 1 

That means the condition <RH.5> is satisfied only for t > 0, i.e. t > 0. In this case, the accelera- 
tion vector reads 

(6.6) X a 



Vcot 2 T + 1 t 

with the radial unit vector u" tangent to the leaf. Condition <RH.6> is just 



(6.7) 



f 3 

r>8 



which means that the condition holds for points which are above the straight line plotted in 



7. Conclusion 



We have shown in this work that the stability of the constraint propagation is heavily influenced 
by the choice of the time slicing, i.e., by the choice of the time coordinate. We derive our res- 
ults by an application of the Routh-Hurwitz criterion to the propagation equation obtained by 
freezing the coefficients in the system of PDEs that describes the propagation of the constraints. 
The system we have treated is known as the Weyl system because it describes the evolution of 
the gravitational degrees of freedom given by the Weyl tensor. We obtain conditions for the 
extrinsic curvature and the acceleration vector of the foliation which have to be satisfied for the 
local constraint propagation to be stable. 

This work can only be considered as a first step towards understanding the various contri- 
butions to the behaviour of the constraints in numerical simulations. For instance, we have 
ignored the effect of the Ricci rotation coefficients to the propagation tensor by focussing ex- 
clusively on the time foliation. 

Furthermore, it is well known that one can add to the evolution equations terms which are ho- 
mogeneous in the constraint quantities. This will generate a subsidiary system whose principal 
part is different from the original one and which will have different stability properties. These 
possibilities have also been completely ignored by us. However, general considerations (see |2J) 
show that such modifications cannot be successful as long as they keep the time reversal sym- 
metry of the Einstein equations. Basically, the argument goes by considering with each initial 
data set the time reversed initial data set. If the perturbations of the constraints are damped in 
one case then by time reversal, they will grow in the other case. Thus, the constraint surface 
cannot be attrative. 

We have also treated the Weyl system in isolation, i.e. decoupled from the rest of the system 
that describes the time evolution of vacuum space-times. In this larger system, there is feedback 
into the Weyl constraints by the constraints coming from the other subsystems and vice versa. 
Again, a more detailed analysis should be made in order to understand these effects. 

Even though we have applied our analysis only to the case of the Weyl system it is clear that 
a similar analysis can also be applied to the standard ADM and the BSSN formulation. Then 
one could compare the analytical findings with the study by Shinkai and Yoneda. It would be 
interesting to see whether the much improved performance of the BSSN system compared to 
the ADM system can be understood from the point of view of this stability analysis. To end, 
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we quickly recapitulate the procedure we have followed for our analysis, and which is to be 
applied to answer the questions given last: 

(1) Starting from the covariant field equations, hyperbolic reduction produces a set of con- 
straint equations and evolution equations. The distribution of the properties of the cov- 
ariant equation into the constraint and evolution equations is decided by the choice of 
gauge which appears as coefficients in the split equations. 

(2) For the constraint quantities whose vanishing indicates fulfillment of the constraint equa- 
tions, propagation equations can be derived using the evolution equations. These are con- 
sidered as equations for the constraint variables propagating on the fixed background 
provided by some solution of the full system of evolution and constraint equations. 

(3) The propagation equations generally couple all of the system's constraint quantities. As 
a first step, feedback of one subsystem into itself can be studied by assuming the other 
constraint quantities to vanish which decouples the subsystem. 

(4) Freezing of coefficients approximates the decoupled system locally by partial differential 
equations with constant coefficients defined on the local tangential space. 

(5) These frozen propagation equations are apt to Fourier analysis yielding propagation equa- 
tions in frequency representation which are just ordinary differential equations with con- 
stant coefficients. 

(6) The location of the complex eigenvalues of the matrix of coefficients (propagation tensor) 
governs the stability behaviour of constraint propagation in this localised picture. Ne- 
cessary for asymptotic stability is that the real parts of the eigenvalues are negative. 

(7) The Routh-Hurwitz criterion is the proper tool to analyse the spectrum with respect to 
stability. It allows to distill algebraic conditions for the gauge parameters under which 
stable constraint propagation is possible. 

These conditions represent necessary and sufficient conditions for locally asymptotic stability 
of constraint propagation within the subsystem under analysis. For the stability of the whole 
system, they are in general necessary conditions. 
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Appendix A. Simplifying assumptions 

In our treatment of the propagation system it has been necessary to make some simplifying 
assumptions in order to bring the equations into a manageable form. Here, we want to dis- 
cuss these assumptions in more detail. Let us start with the covariant form of the propagation 
system <5.H 

(A.l) dfc = ~ tiecrtVJ* + \Uca*X a f h + \Xc h fh ~ jX^c 

This equation is given on the space-time manifold AA which, we imagine, has been foliated into 
space-like leaves If given by t = const, for some global time coordinate t. The time-like vector 
field t a is chosen to be the time-like unit-vector field to the foliation. In addition, we suppose 
that three further space-like unit-vector fields e" (z = 1, 2, 3) have been chosen in order to form a 
complete tetrad field with = t". The space-like members of the tetrad are necessarily tangent 
to the leaves of the foliation. Let (t a = w®,w' a ) be the dual basis. Then we can write every 
tensor field with respect to this basis. Thus we have e.g., 

(A.2) T c = Tiuj c . 

Inserting these expansions into the equations and taking components yields for I — 1, 2, 3 
(A.3) dF, - A k ,n = -i\t{ nn (dmfn - YmK^k) + i^lmnX" 1 ^ 1 + \ X] m T m - \ X T h 

Here, d k := denotes the directional derivative along the tetrad vector e£ while d is the 
directional derivative along t" . Introducing the familiar (3+l)-split we may write 

(A.4) I- = Nd + N k d k 

at 

with the lapse function N and the shift vector N k e k . 

The functions A' c ; and y K are Ricci rotation coefficients with respect to the tetrad defined by 

d4 = A\e>i di4 = ri k el. 

They characterise the behaviour of the spatial tetrad vectors. Geometrically, the functions A k t 
determine how the spatial triad is transported from one leaf of the foliation to the next. In the 
formulations we are interested in, they are considered as gauge source functions, i.e., they can 
be prescribed freely as functions on AA . Here, we will assume that they in fact vanish. This 
amounts to moving the spatial vectors by Fermi- Walker transport along the integral curves of 
t a , the world lines of the Eulerian observers attached to the foliation (i.e., those observers for 
which the current leaf is the manifold of simultaneity). 

Now we fix some event p £ AA which will be the point on which we localise the equation. 
Let fo = t (p) be the value of the time coordinate at p and let Zf be the leaf through p. Let us 
introduce normal coordinates inside Lf g centred at p with respect to the metric h a \, and choose 
the spatial triad at that point to agree with the coordinate vectors. Then only at the point p we 
have Y m k n {p) = 0. After localisation we have dA.3i with all the coefficients being frozen at their 
value at p. This results in 

(A.5) d t F, - N k d k ^, = ~ (-1 £/"'"3 m ^ + 3i ei m nX m F n + x"'^n - 3x?i) ■ 

Admittedly, this procedure of removing the functions Y m k n (p) is somewhat brutal and not quite 
consistent. However, it is the best that we can do if we want to study the isolated influence of 
the time foliation on the stability properties. 
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Appendix B. The generalised Routh-Hurwitz criterion 

The appropriate tool to answer questions of stability of ODEs is the Routh-Hurwitz criterion, 
also known as Bilharz criterion, which originates from stability theory. Presentations of the 
criterion can be found in Gantmacher 1 17J and Parks and Hahn [22 J. Our representation follows 
the lines of 1171 : 

Theorem 1 (Routh-Hurwitz). Let f be a complex polynomial of degree n with 

f(iz) = b z n + b x z n ~ l H +b n + i (a z n + a x z n ~ x + •••+< 

with a,-, b\ real and without loss of generality ^ (otherwise assign f — > if, which does not change 
its roots). Now define the 2p-dimensional Hurwitz determinant 



V 



2p 



a Q a\ 

bo h 

a 

b 



«2p-l 

a 2p-2 
blp-2 



where p = 1, 2, . . . , n and a^ = b^ = Ofor k > n. Further let both of the following real polynomials 



a z n + a x z n ~ x + - 



and z i — ► b§z n + b\Z 



n-1 



a„ 
b„ 



be coprime, which is equivalent to T>m ^ 0. Then the number of roots of the polynomial f, which are 
located in the right complex half-plane {5H(z) > 0}, is given by 



(B.l) 



V(l,V 2 ,V ir ...,V 2n ), 



with V[x\, X2, ■ ■ ■ , x„) denoting the number of changes of sign in the sequence (x\, x 2 , ■ ■ ■ , x„). In case 
some of the determinants in (|B.l^ vanish, then for every section of vanishing determinants of length p 



V 2h ^0,V 2 



(h+1) 



V 



2{h+p) 



0,V 2(h+f , +1) ^0 (h, p£N) 



one has to set: 



(B.2) V(V 2h , T>2(h+l)' "D 2 (h+py ^2(h+p+\)) 



P+i 

? + 2~ £ with e 



p. 



(-1)4 sign ^ 



2(ft+p+l) 



2h 



ifp odd 
ifp even 



The proof is rather extensive and can be found in the mentioned literature, the basic idea how- 
ever will be outlined in the following section. From the theorem follows, that equivalent to 
asymptotic stability is, that all Hurwitz determinants are strictly positive. 

B.l. Mathematical background. The basic idea behind the remarkable Routh-Hurwitz criterion 
is the argument principle: Every polynomial / of degree n has exactly n complex roots oq and can 
therefore be written as a product of its elementary divisors: 



(B.3) 



/( z ) = e a ' 2 ' = a " ri( z ~ a i) 

i=0 i=l 
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Then the argument of / constitutes additively from the contributions of the several elementary 
divisors: 

(B.4) arg/(z) = arg (^a n JJ(z - = £ arg(z - a t ) + arga„ 

Now choose a closed, non-self-intersecting path ^ in the complex plane and track the change 
of arg/(z) = E? = i ar g( z ~~ <*i) + ar g fl n while travelling once around c € in positive orientation 
(s- fig-EJ. 




FIGURE 4. The argument principle 



As the change of argument sums up from the contributions of each root, one can look at each 
root individually: 

If the considered root lies outside of the domain enclosed by (e.g. CX2 in fig. 0), then the 
argument of (z — oci) first grows a little, then decreases and finally increases so that the net 
difference exactly vanishes. Roots out of ^ therefore contribute nothing to the change of argu- 
ment. 

On the contrary, if the root under consideration lies inside c € (as ot\ of fig. 0| does), then the 
argument of (z — <X\) grows continually, picking up an increase of 2zr for one revolution. 

So the growth of argument counts the number of roots I inside ^ , counting multiple roots 
according to their multiplicity: 

(B.5) A^arg/ = 27r/ 

This connection between change of argument and the location of roots can now be employed 
to count the roots in a complex half-plane. Let I and r be the number of roots in the left and 
right half-plane, and n = I + r, i.e. no roots lie on the imaginary axis. Now one can construct a 
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path <# :-- 
R —> oo. 
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fl^2 out of two segments (s. fig-EJ, which encloses the left complex half -plane for 




FIGURE 5. A path enclosing the left complex half -plane for R — ► oo 



Calculating the change of argument for the individual segments, the half -circle ^2 ror R ~^ 00 
always counts all roots, regardless of their location: 



(B.6) lim A V2 arg/ = Km £ A^ arg(z - aci) = £ Jim arg(R(^ - f )) 



nn 



R^oo 



TT 

' 2 



The complete path c if however counts only the roots in the left complex half-plane according to 
the argument principle: 



(B.7) 



lim A^ arg/ = In I 

R^oo 



Thus the contribution of ^ along the imaginary axis is given as the difference of these two 
amounts: 

(B.8) lim A'g ar g/ = li m &v ar g/ — lim Acg arg/ = n (21 — n) = n (n — 2r) 

R—>oo R—>oo R— >oo 



lim Atf arg / = Aarg/(iz) 

R— >oo 



This can be written as 
(B.9) 

Representing the polynomial / by 

(B.10) /(iz)=:/jR(z)+// 3 (z) 



z^oo 

z— > — OO ' 
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with real polynomials fa and fo, and using 



(B.ll) arg/(t'z) = arctan 



M(z) ' 



now allows to state the Leonhard-Michailov criterion: 

All roots are located in the left complex half-plane (r = 0), if and only if 

(B.12) arctan -^ Z ' 



increases exactly by n n when travelling along the real axis from z = — oo to z 
That means, one has to track the change of angle of the vector 



(B.13) 



Mz) 



when walking from z = — oo to z = +00, and count how often the vector circles around the 
origin. This is a rather cumbersome method. Fortunately, it is possible to find an algebraic 
version of this criterion. The following shall outline the necessary procedure. For a proof, we 
refer to literature, e.g. |22 section 1.2, p. 10ff]. 

As a first step, the growth of arctan by multiples of n is directly related to the number 

f (z) 

of jumps between —00 and +00 of j^y- For every increase by re, the fraction jumps exactly 
once from +00 to —00. 

A magnitude which counts such jumps is the Cauchy index I%g which amounts to the number 
of jumps of the function g from —00 to +00 minus the number of jumps from +00 to —00, when 
tracing g from a to b. 

Considering rational functions the Cauchy index can be calculated from the Sturm se- 
quence if s\ and S2 are coprime and the degree of Sj is greater than that of S2. The Sturm sequence 
is generated by the Euclidean algorithm: For two elements s t -_i and s; of a Sturm sequence, the 
next element s,- +1 is the negative residual of the polynomial division (s ; _i : s,), i.e. 



(B.14) s,_i = qi_ x si 



n+i, 



where the quotient polynomials (<ft) are of no further interest. As the degree is decreasing from 
element to element, the Sturm sequence terminates with a constant polynomial s„, which is non- 
vanishing if and only if s-[ and S2 are coprime. In this case furthermore the Sturm theorem for the 
Cauchy index of the rational function holds: 

(B.15) I*A = V(a)-V(b), 

with V(z) denoting the number of changes of sign in the sequence (si(z),S2(z), ■ ■ ■ ,s m (z)). 
Thus one evaluates the elements of the Sturm sequence at the position z and then counts the 
number of changes of sign. 

For the Leonhard-Michailov criterion one is interested in the Cauchy index therefore 
the Sturm sequence must be evaluated at z — > —00 and z — > +00. For these limiting values 
only the highest order coefficients in the elements of the Sturm sequence play a role. They can 
be expressed in terms of special determinants, which are composed out of the coefficients of fa 
and fa. 

Finally this leads to the Routh-Hurzvitz criterion given above, which therefore represents an 
algebraic version of the Leonhard-Michailov criterion. 
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